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Abstract. We show how to ehminate the error caused by an incorrectly modeled bound- 
ary in electrical impedance tomography (EIT). In practical measurements, one usually lacks 
the exact knowledge of the boundary. Because of this the numerical reconstruction from 
the measured EIT-data is done using a model domain that represents the best guess for 
the true domain. However, it has been noticed that the inaccurate model of the bound- 
ary causes severe errors for the reconstructions. We introduce a new algorithm to find a 
deformed image of the original isotropic conductivity based on the theory of Teichmiiller 
spaces and implement it numerically. 
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1. Introduction. We consider the electrical impedance tomography (EIT) problem, 
i.e. the determination of an unknown conductivity distribution inside a domain, for ex- 
ample the human thorax, from voltage and current measurements made on the boundary. 
Mathematically this is formulated as follows: Let f2 C be the measurement domain, and 
denote by 7 = (7*-^) the symmetric matrix describing the conductivity in Q. We assume 
that the matrix has components in L^{Q) and that it is strictly positive definite, that is, 
for some c > we have {^,^{x)0 > c||^|p for all x E fl. The electrical potential u satisfies 
in n the equation 

(1.1) V-7Vm = 0. 

To uniquely fix the solution u it is enough to give its value on the boundary. Let this be 
u\dn = / £ H^/'^{dQ) where H^/'^{dfl) is the Sobolev space. Then (jl.ip has a unique weak 
solution u G H^{Q). 

Our boundary data is the map that takes the voltage distribution / on the boundary 
for all / to the corresponding current flux through the boundary, u ■ 7VM, where u is the 
exterior unit normal to Q. Mathematically this amounts to the knowledge of the Dirichlet- 
Neumann map A corresponding to 7, i.e. the map taking the Dirichlet boundary values to 
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the corresponding Neumann boundary values of the solution to 



dxj 



an 



This defines a bounded operator : H^^'^^dVt) H ^/^((9f2). The symmetric quadratic 
form corresponding to A^, 

(1.2) A^[h,h]:= / hA^hdS = / Vu--fVudx 

JdQ JdQ 

equals in physical terms the power needed to maintain the potential h on dQ. 

When 7 is a scalar valued function times identity matrix, we say that the conductivity 
is isotropic. As usual, conductivities that may be matrix-valued are referred as anisotropic 
conductivities. The EIT-problem is to reconstruct 7 from A^. The problem was originally 
proposed by Calderon j6J and then solved in dimensions three and higher for isotropic 
smooth conductivities in ^H]- The two dimensional case that is relevant to us was solved 
by Nachman ^3] for isotropic conductivities assuming 7 G W'^'^, p > 1 and then finally 
for general L°*^-smooth isotropic conductivities by Astala and Paivarinta in a celebrated 
paper jl]. 

The conductivity equation is invariant under deformations of the domain Q in the fol- 
lowing sense. If F is a diffeomorphism taking Q to some other domain Q, then u o will 
satisfy the conductivity equation in Q with conductivity 



detF'(y)| 



y=F-^x) 



where F' is the Jacobi matrix of map F and u is a solution of V- 7VU = in ^2. We say 
that 7 is the push forward of 7 by F and denote it by 7 = ^^,7. Note that all this is well 
defined for general matrix valued 7. For us the starting point is the trivial observation that 
even if 7 is isotropic, the deformed conductivity 7 will not in general be isotropic. The 
boundary measurements are invariant: When / : dQ — dQ is the restriction of F : Q ^ Q, 
we say that A = f^A-y, 

{{f.A,)h){x) = {A,{h o /))(y)|^^^_,(,.) , h G H'/\dQ) 

is the push forward of A^ in /. As seen in [T7], it turns out that /*A^ = Ap^^. 

The facts that the anisotropic conductivity equation and the boundary measurements 
are invariant has the important consequence that the EIT-problem with an anisotropic 
conductivity is not uniquely solvable, even though the isotropic problem is, see jT7j . 

In practice when solving the EIT-problem in a given domain Q, one typically seeks for 
the isotropic conductivity that minimizes 



'1-4) IIA^ea. -A^|p + a||7||^ 



for 7 defined in terms of some triangulation of Q as e.g. a piecewise constant function and 
II ■ ||x is some regularization norm jTU]. Here A^eas is the measurement of the Dirichlet- 
Neumann map that contains measurement errors. 
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In practice, one of the key difficulties in solving the EIT problem is that the domain Q 
may not be known accurately. It has been noticed that the use of slightly incorrect model for 
fl, i.e., a slightly incorrect model of the boundary causes serious errors in reconstructions, 
see e.g. [HICIIH]- As an example, consider the EIT measurements of pulmonary function 
from the human thorax. The measurement electrodes are attached on the skin of the 
patient around the thorax. In principle, an exact parameterization for the shape of the 
thorax could be obtained from other medical imaging modalities such as magnetic resonance 
imaging (MRI) or computerized tomography (CT). However, in most cases these data is 
not available, and one has to resort to some approximate thorax model. Further, the shape 
of the thorax varies between breathing states, and it is also dependent on the orientation 
of the patient. Thus, the thorax geometry is known inaccurately even in the best case 
scenarios. 

In this paper our aim is to propose a method to overcome the problem that the boundary 
and its parameterization are not exactly known. The set-up of the problem we consider is 
the following: 

We want to recover the unknown conductivity 7 in i7 from the measurements of Dirichlet- 
to-Neumann map, and we assume a'priori that 7 is isotropic. We assume dfl and are 
not known. Instead, let Q^, called the model domain, be our best guess for the domain 
and let fm '■ dQ — >■ dflm be a diffeomorphism modeling the approximate knowledge of the 
boundary. As the data for the inverse problem, we assume that we are given the boundary 
of the model domain, dflm and the map A^ := (/m)*A^ on dQm- Note that we have 
simplified the problem by assuming that the only error in A^ comes from the imperfect 
knowledge of the boundary. 

This set-up is motivated by the fact that the quadratic form corresponding to A^, 

Am[Kh]= [ hAmhdS= [ {hofJA^{hofJdS, heH^'\dn^) 

equals physically to the power needed to maintain the potential ho on dQ. 

Since A^ usually does not correspond to any isotropic conductivity because of the de- 
formation done when going from the original domain Q to Q^, we obtain an erroneous 
solution 7 when solving the minimization problem (|1.4|) . This means that a systematic 
error in domain model causes a systematic error to the reconstruction. In particular, local 
changes of the conductivity often give raise to non-localized changes in reconstructions due 
to the above systematic error. Thus the spatial resolution of details of reconstructions are 
often weak. This is clearly seen in practical measurements, see e.g. jH|. 

We note that one could forget in solving of the minimization problem (jl.4j) the assump- 
tion that 7 is isotropic, and find the minimizer in the set of anisotropic conductivities. 
However, the anisotropic inverse problem has non-unique solution, and as the minimiza- 
tion problem is highly non-convex, the minimization would be hard and, as usual, forgetting 
existing a'priori information makes the solution significantly worse. 

To formulate our main results, let us define certain concepts. We start with the maximal 
anisotropy of an anisotropic conductivity. 

Definition 1.1. Let -^^^[x) he an L°°{Q)- smooth matrix valued conductivity in Q and let 
Xi{x) and X2{x), Xi{x) > X2{x) be the eigenvalues of matrix '-j^^{x). We define the maximal 
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anisotropy of a conductivity to be K{'j) given by 



K (•y) = sup K {•y , x) , where K{y^x) 



Ai(x) 



We call the function K{'-j^x) the anisotropy of y at x. Here sup denotes the essential 
supremum. 

Sometimes, to indicate the domain Q, we denote -^'(7) = -^"^(7). As a particularly 
important example needed later, let us consider the conductivity matrices of the form 



where A > 1 is a constant, r]{x) G M+ is a real valued function, Re^x) is a rotation matrix 
corresponding to angle 9{x), where 



We denote such conductivities by 7 = yx,e,r)- These conductivities have the anisotropy 
K{'~f,x) = c\ = {X^^"^ — 1)/(A"^/^ + 1) at every point and thus their maximal anisotropy is 
K = cx. We call such conductivities 7 uniformly anisotropic conductivities. 
The following theorem is the main result of the paper. 

Theorem 1.2. Let Q be a bounded, simply connected C^'" -domain with a > 0. Assume 
that 7 G C°'"(fi) is an isotropic conductivity and its Dirichlet-Neumann map. Let 
be a model of the domain ( which is assumed to satisfy the same regularity assumptions as 
fl), and fm '■ dVt dQ^n be a C^'°'-smooth diffeomorphism. 
Assume that we are given dQm o-nd Km = (/m)*A^. Then 

1. There are unique X>1, 6 E L°°{Qjn, S^) and rj G L°°{Qrn, ^+) such that the conduc- 
tivity 7 = 7A,6»,r? satisfies = A^. 

2. If ^2 is an anisotropic conductivity in Qm such that A^^ = Am then -^'(72) > -^^(7)- 
Moreover, -^'(72) = -^^(7) if and only if 72 = 7. 

Theorem 11.21 can be interpreted by saying that we can find a unique conductivity in Qm 
that is as close as possible to being isotropic. 

The proof of Theorem 11.21 is based on the theory of quasiconformal maps. There are 
several equivalent definitions for these maps, and we will present the one based on a partial 
differential equation (Beltrami equation) in Section 2. However, the quasiconformal maps 
have also a geometric definition. Indeed, they are generalizations of conformal maps that 
take infinitesimal disks at z to infinitesimal disks at f{z), and the radii gets dilated by 
\f'{z)\. Analogously, a homeomorphic map is quasiconformal on a domain Q if infinitesimal 
disks at any z E Q get mapped to infinitesimal ellipsoids at f{z). The ratio of the larger 
semiaxis to smaller semiaxis is called the dilation of / at z, and the supremum of dilatations 
over Q is the maximal dilation. This dilatation of infinitesimal discs is in fact the reason why 



(1.5) 
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isotropic conductivities change to anisotropic ones in push forwards with quasiconformal 
maps. 

The crucial fact that we use proving Theorem II .21 is a result of Strebel [T^, that roughly 
speaking says that among all quasiconformal self-maps of the unit disk to itself with a given 
sufficiently smooth boundary value there is a unique one with the minimal maximal dilation. 
This will yield that corresponding to the given boundary modeling map : dQ — * dQ^n 
there is a unique map F : Q ^ Qrn having the minimal maximal dilation. We will show 
that this leads to the following result: 

Proposition 1.3. Let Q, 7, Qm, (^nd fm satisfy assumptions of Theorem M.'A Then there 
is a unique map F : Q ^ Qm, depending only on f : dQ —>■ dQm such that for the uniformly 
anisotropic conductivity 7 corresponding to 7 in Theorem \l.^ we have 

(1.6) det(7(x))i/2 ^ ^(F-i(a;)). 

Proposition II .31 can interpreted as saying that solving the minimization problem p.4|l with 
a = in the class of conductivities ^x,e,n we can find the function (det7(x))^/^ in f2m 
that represents a deformed image of original conductivity 7 in the unknown domain Q and 
the deformation depends only on the error made in modeling the boundary, not on the 
conductivity in Q. 

In particular, this turns out to be useful as local perturbations of conductivity remain 
local in reconstruction: if we consider one fixed boundary modeling map fm '■ dQ dQrn 
but two isotropic conductivities 71 and 72 = 71 + a in Q, then the reconstructions 71 and 
72 obtained by Theorem 11.21 corresponding to 71 and 72 satisfy 

det(72(a;))'/' - det(7i(a;))^/' = (r{F'\x)). 

Remark 1. Theorem 11.21 holds for anisotropic conductivities in the sense that for each 
C^'"(f2)-smooth anisotropic conductivity 7 there is a unique conductivity 7 = 7A,e,J7 such 
that (/m)*A^ = A;^. However, for anisotropic conductivities Proposition 11.31 does not hold 
as the map F depends on the conformal class of 7. 

The paper is organized as follows. In Section 2 we consider isotropization of anisotropic 
conductivity using a diffeomorphism and pay close attention to the smoothness required 
from 7 and Q and introduce the necessary background from the theory anisotropic inverse 
problems. We apply this in Section 3 in proving main results using the existence of a 
Teichmiiller mapping. In Section 4 we consider physically realistic measurements, i.e., so- 
called complete electrode model. The numerical implementation for the complete electrode 
model is then described in the last sections. 

2. Quasiconformal maps and solvability of inverse problem with anisotropic 
conductivity. It is a classical result that every Riemannian surface is locally conformal 
to a Euclidean plane: this corresponds to choosing the coordinate system to be isother- 
mal. Similarly, every anisotropic conductivity matrix can be transformed to an isotropic 
conductivity. We identify the plane with complex plane C. 
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Lemma 2.1. Let Q be a bounded, simply connected C^'"' -domain with a > 0. Assume 
that 7 G C^'°'{Q) is an isotropic conductivity. Then there is a C^'"' -smooth diffeomorphism 

F -M^n, n = F{n) C C such that 

(2.7) F,7 = 

where F^'y is defined by il.!^} . and (3 is the identity matrix multiplied by a C^'"' -smooth 
scalar function. Moreover, 

(3= (det7oF-i)i/2j_ 

The proof of this result is well known, but as smoothness of F is crucial later, we give 
the proof for the convenience of the reader. 

Proof. The equation ()2.7|) a'priori a nonlinear system for the derivatives of F. However, 
in two dimensions this equation completely linearizes, and is equivalent to the Beltrami- 
equation 

(2.8) dF = ndF, 

where the complex derivatives are d = |(^ — '^^); d = |(^ + i-§^) and the Beltrami 
coefficient /i = fipiz), called also complex dilatation is given by 

n^ + ^22 - 2i7i2 

(2.9) /X 



7ll + 722 + 2a/7ii722 - 7l2 
The function /i has the crucial property that it is strictly less than one in modulus: 

(2.10) sup |/i(2)| < 1. 

Let us extend the conductivity matrix 7 (a'priori only defined in Q) to the whole plane to 
be the identity matrix outside fl. Similarly, fi is extended outside Q by zero. 

Next we consider how to solve the Beltrami equation, and for this we consider it in the 
whole plane. In order to have a unique solution we fix the behaviour of F at infinity. Thus, 
consider 

(2.11) dF{z) = ij{z)dF{z) inC, 

F{z) = z + h{z), 

lim h{z) = 

where /i is a compactly supported L°°-function satisfying ()2.1()|1 . This problem has unique 
solution F & when p is close enough to 2 and —2/p<6<l — l/p. For the proof of this 
see, for example |2] or [TTj . The proof is based on the fact that ()2.1H) can be written as an 
integral equation 

(2.12) ^^^)-2^X^3^^"^^) = " 

where da{C,) is Euclidean area in C (or R^). As ||yu||oo < 1, it turns out that the left 
hand side of equation ()2.12j) is of the form of the identity plus a contractive operator in 
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Sobolev space W^'^^fl), with appropriate p, and thus equation (j2.12|) can be solved by an 
apphcation of the Neumann-series argument. 

Using interior Schauder estimates for equation ()2.1H) . we see that if 7 and thus /i are 
C'^'^'-smooth, the solution F has to be locally C^'"-smooth in C, in particular in Q. Using 
formula ()1.3p we see that F^,7 is C'^'^'-smooth in closure of f2. □. 

In general, any solution F : Q ^ Q to Beltrami equation for fi satisfying (j2.1(Jj) and 
for which F G H^{Q) is called quasiregular. If a quasiregular map F : f2 — > f2 is homeo- 
morphism, it is said to be quasiconformal. The quasiconformality can be defined also in 
geometrical terms, see [H]. 

Next we recall the recent results for inverse problems for anisotropic conductivities 7. 
Let us consider a class of conductivities in Q, given by 

S(7) = {F*7 \ F -.n^nis homeomorphism , F, F"^ G H^{n; C), = /}, 

that is, ^(7) is the equivalence class of the conductivity 7 in push forwards with boundary 
preserving diffeomorphisms. Then A^- = for all a G ^(7). By 0, the converse is true, 
that is, if 0" is a strictly positive definite L°°-conductivity and A^- = A^, then a G ^(7). In 
other words, A^ determines the equivalence class ^(7). Note that diffeomorphism F : Q 
fl such that F G H^{Q; C) and F\gQ = I is quasiconformal. 

3. Proof of main results We start with the proof of Theorem 11.21 Note the fact that 
the conductivity 7 in i7 is isotropic will not be used in the proof at all. First we show that 
we can assume that Qm is the unit disc D C C. 

To prove this, let fm '■ dfl dQm be the boundary modeling map and 7 a C^'"-smooth 
(also possibly anisotropic) conductivity in Q. 

Our first observation is that as flm is a simply connected domain, it follows from Riemann 
mapping theorem that it can be mapped to unit disc © conformally. Moreover, as flm is 
C^'"-smooth domain it follows from Kellog-Warschawski theorem ^1 Thm. 3.6], that the 
Riemann-map can be chosen to be be a C^'"-diffeomorphism Fq : flrn ^ such that 
Fq : Qm ^ D is conformal. Thus, if a is some C°'"-smooth anisotropic conductivity in Q^n, 
we have that (Tq = (Fo)*o" is C°'"(ro)-smooth conductivity. 

Second, we observe that the uniformly anisotropic conductivity 7A,e,»7 of the form ()1.5|) 
in Qrn changes under (Fq)* to a uniformly anisotropic conductivity (Fo)^,7a,6i,j7 = 7a,6»o,»?o i^ 
D such that rjo = r] o F^^. 

Third, we see that as Fq : Qm — > © is conformal, the maximal anisotropy of (Fo)*cr and 
cr satisfy 

Ko{{Fo),a)=Kn{a), 

that is, the maximal anisotropy is preserved in conformal transformations for any a. 

Fourth, if /q = Fq\qq^ then A^^ = (/o)*Ao-. Also, we see that the our data is invariant 
in the change of the model domain in the sense that (/m)*A^ = (/o)*((/m)*A^) where 
fm = foo fm ■ dil ^ dB. 

These four observations yield that it is enough to prove the assertion in the case when 
Qm = Indeed, changing Q^n to D with Fq keeps the boundary measurements, the 
smoothness of objects, the maximal anisotropy as well as class of uniformly anisotropic 
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conductivities invariant. More precisely, we can replace the boundary modeling map fm by 
the map fm = fo°fm- 

Thus, let us return proving Theorem 11.21 in the case when Qm = D- Let fm '■ dQ dV> 
be the boundary modeling map that is a C^'"-smooth diffeomorphism and 7 G be 
an isotropic conductivity with Dirichlet-Neumann map A^. 




Let Fm be some C^'"(f2)-diffeomorphism F^ : — D such that Fm\dn = fm- There are 
many ways to construct such map, and for convenience of the reader, we present one simple 
way. Let G : i7 — * D be a Riemann-map. By Thm. 3.7], G has C^'^'-extension G : Q 
D. Let = fmoG"^ -.OB^dB and = \z\ exp{i{\z\^aTg{(f){z/\z\)) + {l-\zf)aTg{z/\z\))) 
be a C^'"-diffeomorphism © ^ D satisfying $|aD = 0. Then Fm can be chosen to be the 
map $ o G. 

Let 7o = (-Fm)*7 be an anisotropic conductivity in D. By Lemma l2.ll there is a C^'"- 
diffeomorphism Fi : Q —>■ Qi such that the conductivity 71 = (-Fi)*7o is C°'"-smooth 
isotropic conductivity in the closure of the C^'"-smooth domain Qi. 

As n 1 is a simply connected C^'"-smooth domain, by Kellog-Warschawski theorem 
cited above there is a conformal map F2 : Qi ^ B> such that F2 : fii — D is a C^'°- 
diffeomorphism. Let ^3 = ^20^1:©-^© and /s = F^Iq^. Note that (-F3)*7o is isotropic 
conductivity in D as F2 is conformal. 

The boundary values of quasiconformal maps D — © are characterized as being the 
quasi-symmetric maps, that is, homeomorphic maps / : d3 dD such that 6{u) = 
arg/(e™) satisfy 

(3.13) <^!rl^^4!rT <k, ioi all u,v e dB, 

9{v) — 9[u) 

with some k > 0, see |llj . 

Let us consider next the map f^ = f^^ : dlD) dlD). Since /s and /3~^ are C^'"-smooth, 
we see that /4 satisfies 

0(u + t) - 9(u) , , 

(3.14) lim -) { = 1 uniformly in m G 9D, 

^ ^ i-o e{u -t)- e{u) 
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and is in particular quasi-symmetric. Thus /4 is boundary value of at least one quasi- 
conformal map. What is more, since f\ satisfies condition (j3.14p . it follows from results 
Strebel |16J, that among all quasiconformal maps having f\ as a boundary value there is 
a unique extremal map Fe in the sense that the L°°-norm of the complex dilatation is 
minimal. More precisely, if F : D — D is quasiconformal map such that Flsn = f^, then 
its Beltrami coefficient satisfies | l/i^^l |j;^oo > ||/iFel|L°° and the equality holds only if F = F^. 
Furthermore, the extremal F^ is a Teichmiiller mapping, i.e., its complex dilatation fip^ is 
of the form 



(p[z] 

(3.15) f^FA^) = Wf^Fjloo T-r-^ 
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where : D ^ C is holomorphic in D, and has thus discrete set of zeros. Note that as 
Ff, need not to be even Lipschitz smooth near zeros of (p. Reader should also note that 
certain assumption on the regularity of the boundary value fm is necessary for existence of 
extremal maps. This will be discussed after finishing the proof. 

Let us now consider how a quasiconformal map F : D — > D with complex dilatation fip 
change maximal anisotropy of conductivities. When a is an isotropic conductivity in D, 
that is, K{a) = 0, one sees that for the anisotropic conductivity a = F^a we have 

K{x, a) = fiF{F'\x)), for x e 1, 

and hence the maximal anisotropy satisfies Kla) = ||/iF||L°°- 

Let now 73 = (^3)^,79 be an isotropic conductivity in D and let de = (Fe)^,73 be an 
anisotropic conductivity in D. Here, F^o F3 o fmlan = fm- In particular, the above shows 
that 

(/m)*A^ = (/4 O /3 O fm)*A^ = {fi O /3)*A^(, = (/4)*A73 = A^^. 

In particular, this implies that inverse problem of finding conductivities cr in D such that 
(/m)*A-y = Aq- has a solution cr = a^- By Section |21 the knowledge of the boundary 
dflm = dB> and the map (/m)*A-^ determines the class S(cre) of conductivities in D. Now 
we can write the class S((Te) also as 

S((Te) = {F,73 : F : D ^ D is homeomorphism, F, F"^ G H\Q; C), F|ao = f^}. 

Since 

K{F^J3) = ||/iF||L«=(D), 

we see that the conductivity cTg = (Fe)*73 corresponding to the extremal map Fg is the 
unique conductivity a in the class S(o"e) that has the smallest possible value of K{a). 

Finally, since \fJ^Fe{z)\ = cq is constant function of 2; e D, and cxe = (Fe)*73 with isotropic 
73, we see that the ratio of the eigenvalues of the conductivity matrix ae{z) is constant for 
z E 3. Thus cTe has the form cTg = 7A,e,?j with cq = (1 — A)/(l + A), = 73 o (Fg)""^, and 
some 6. This proves Theorem 11.21 □ 

As noted above, the construction of 7 in the above proof did not use the fact that 7 is 
isotropic at all. This shows Remark 1. 
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Next we proof Proposition 11.31 

Proof. Consider isotropic conductivities 71 and 72 in Q. In sequel, we use notation of the 
proof of Theorem ll.2[ By definition, fm determines a map Fm- The construction of the map 
Fi is based on the Beltrami coefficient of the conductivity. Clearly, the Beltrami coefficients 
for the conductivities (i^m)*7i and (-Fm)*72 coincide, and thus Fi and fli can be taken to 
be the same for both 71 and 72. The maps F2, F3 and Fe are constructed by using dQi and 
Fi, and thus they coincide for 71 and 72. Since in general det(F*7)(x) = det(7(F~-'^(x))), 
this proves Proposition 11.31 □ 

As noted above, certain assumption on the regularity of the boundary values are neces- 
sary, for there are counterexamples to the uniqueness, for example the so-called Strebel's 
chimney. The current state of the uniqueness question can be found in [H]. 

We note also that if in formula (j3.15p the function has zeros in fl, then fip^ has a 
singularity of type {z — /{z — zqY , and this could affect the behaviour of the recon- 
struction algorithm we propose in a way to explained later. However, in all the numerical 
examples we have tested these difficulties do not appear, probably since our deformations 
are relatively small. 



4. Electrode model In the numerical simulations below we have used so called com- 
plete electrode model which is a certain finite dimensional approximation of Dirichlet- 
to-Neumann map. This model is chosen as it is an accurate model for the measurements 
made in practice. As noted before, in experimental measurements one places the mea- 
surement electrodes on the boundary, e.g., the skin of the patient, without knowing exact 
parameterization of the boundary. Thus this model is a paradigm of the case when the 
boundary is unknown. 

To define the electrode model, let ej C dQ, j = 1, . . . , J be disjoint open paths modelling 
the electrodes that are used for the measurements. Let u solve the equation 

(4.16) V-7Vt; = in fi, 

(4.17) Zjiy-fVv + v\e^ =Vj, 

(4.18) ^-TVt^lanvu/^^e, =0, 

where Vj are constants representing electric potentials on electrode ej. This models the 
case where electrodes ej having potentials Vj are attached to the boundary, Zj is the contact 
impedance between electrode Cj and the body surface, and the normal current outside the 
electrodes vanish. By (|4. 16114. 1'Hj) has a solution u G H^{Q). The measurements in 
this model are the currents observed on the electrodes, given by 

1 f 

Ij = - — j- / z/- 7Vt'(x) ds{x), j = 1, . . . , J. 

Thus the electrode measurements are given by map E : M"' —>■ M"^, E{Vi,...,Vj) = 
(Ji, . . . , Ij). We say that E is the electrode measurement matrix for {dQ, 7, ei, . . . , ej, 
Zi, . . . ,zj). Let Q and Q be C^'^'-smooth domains. We say that / : dQ dQ is length 
preserving on Uj^^Cj if ||D/(r)|| = 1 for x G ^j=iej where r is the unit tangent vector of 
dQ. 
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Proposition 4.1. Let and be C^^"" -smooth domains and F : Vt ^ Vt he a C^'°- 
diffeomorphism, Cj C dQ be disjoint open sets, and 'y be a conductivity on Q. Let f = 
F\dn, = /(cj) C dfl and 7 = (-F)*7. Assume that f is length preserving on Uj^^Cj. 
Then the electrode measurement matrices E for {dQ, 7, ei, . . . , cj, zi, . . . , zj) and E for 
{dQ, 7, ei, . . . , ej, zi, . . . , zj) coincide. 

Proof. We start with an invariant formulation of electrode measurements E. For this, 
let R be the Robin-to- Neumann map given by Rf = u- 7VM|an where u is solution of 

(4.19) V-7Vm = inn 

zu-'yVv + rjvldn = h, 

where z = z{x) is C°°{dQ) fuction such that z\ej = zj and 77 = J2j=i Xe^ (x), where Xej is the 
characteristic function of electrode Cj. Note that if the boundary and the contact impedance 
are known, the Robin-to-Neumann and the Dirichlet-to-Neumann maps determine each 
other, that is, they represent equivalent information. 

Consider now the bilinear form corresponding to liner maps E : M"' x M.'^ M and 
R : H-^/^{dn) X H-^/\dn) R given by 

E[V,V] = y2iEV)jVj\ej\, R[h,h]= [ {Rh)hds. 

j=l Jd^l 

Let S = span(xe, : J = 1, . . . , J) C H-^l\dn) and define M:V = {V,)j^^ ^ ^i^-. 
to be a map M : —>■ S. Then 

(4.20) E[V, V] = R[MV, MV]. 
Moreover, for h = MV with some V G M"', we have 

{A.21) R[h,h] = / {u + ziy-'y\/u)iy-'y\/uds = / 'yWu-Wudx + / z |z/- 7Vm|^ rfs 
Jon Jn Jan 

where u solves ()4.19|) . Above the integral over Q is invariant in coordinate deformations. 
Note that in the above formula the integral over the boundary is not coordinate invariant. 

Let E be electrode measurement matrix for 7 in f2 with electrodes ej = f{ej) and let R be 
the Robin-to-Neumann map for 7 defined analogously to (|4.19p . Since / is length preserving 
on the electrodes, we see using (|1.3p that z/-7Vm(x) = z/- 7Vm(/(x)), for m = m o and 
X G dfl, and thus we can see from 1)4.211) that 

R[h,h] = R[hof-\hof-\ 

for h,hG H^^/'^{dVt) supported in the closure of IJj=i^i- Thus for the map M : V 

Y^'^iXZj we have by formula dOIlD that E[V,V] = R[MV,MV]. Combining this and 
(|On|l we obtain 

E[V,V'] = E[V,V']. 
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In particular, this implies that the matrices E and E coincide. □ 

In particular, it the case where Q is the model domain and f = fm '■ dQ dQm is 
the model map for the boundary, the assumption that / is length preserving on electrodes 
means the very natural assumption that in electrode measurements the paramitrization 
of electrodes are known. Then by Proposition 14.11 the electrode model discretization E 
of equals the corresponding discretization E of (/m)*A^. Summarizing, the electrode 
measurements does not change if we have modeled the geometry of the boundary incorrectly 
but the electrodes are modeled correctly. 

5. Numerical examples The performance of the proposed method is evaluated by test 
cases with simulated EIT data. First, in Section f5.ll we briefly discuss the discretization 
and the computational methods that are used, and the results are then given in Section 

5.1. Discretization and notation The numerical solution of the forward model is based 
on the finite element method (FEM). The variational formulation and the finite element 
discretization of the electrode model ()4.16H4.l'5j) in the case of isotropic conductivities have 
been previously discussed e.q. in ^U]. The extension of the FEM-model to the anisotropic 
case is straightforward, the details will be given in a subsequent publication. 

For the functions r]{x) and 6{x) in equation (jl.5j) we use piecewise constant approxima- 
tions that are defined on a lattice of regular pixels. Thus, we have 

M M 

(5.22) ^ = J2r],x^{x), e = Y,Ga'^{x) 

i=l i=l 

where Xi is the characteristic function of the i^^ pixel in the lattice. Within the discretization 
()5.22|) . the parameters r] and 9 are identified with the coefficient vectors 

r/ = (r/i,r/2,...,r/M)''GM^ 
= {61,62, ... , 6m)^ G 

and A is a scalar parameter. Note that as 7A,r?,e = lx',Ti,e', where A' = 1/A and 6'{x) = 
6{x) + 7t/2, we can assume in looking the minimizing uniformly anisotropic conductivity 
that A gets values A > 0. 

In practical EIT devices, the measurements are made such that known currents are 
injected into the domain Q through some of the the electrodes at dfl, and the corresponding 
voltages needeed to maintain these currents are measured on some of electrodes. Often, 
voltages are measured only on those electrodes that are not used to inject current. Thus, 
measurements made give only partial information on the matrix E. To take this in to 
account, we introduce the following notation for the discretized problem. We assume 
that the EIT experiment consists of a set of K partial voltage measurements, V^^\ j = 
1, . . . , K. For each measurement, consider a current pattern j = 1, . . . , K such that 
J2i=i ^e^ — 0- Typically, the corresponding measurements are the voltages (potential 
differences) between pairs of neighboring electrodes. Let us assume that the measurement 
vector V^^^ corresponding to the current pattern I^^^ consist of L voltage measurements. 
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I.e., we have G M^. Thus, we write V^^^ = PjE ^I^^'' + e^^\ where E is the electrode 
measurement matrix, random vector e^^^ models the observation errors and Pj : M"^ M.^ 
is a measurement operator that maps the electrode potentials to measured voltages. 

In the inverse problem, the voltage measurements V^^\ V^'^\ . . . , V^^^ are concatenated 
into a single vector 

V = {V^'\ V^^\ . . . , V^^'Y G M^, N = KL. 

For the finite element based discretization of the forward problem U : t— > M^, we 
use the notation 

t/(r/, A) = (f/«(r/, A), U^^\^, ^, A), . . . , [/^(r/, ^, A))^ G M^, 

respectively. Here, U^^^{ri,9, X) = PjE~'^{r],9, X)I^^^ G M"^ corresponds to partial voltage 
measurement with current pattern I^^^ and conductivity 7r),e,A- 

Using the above notations, we write the discretized and regularized version of our inverse 
problem as finding minimizer of the functional 

(5.23) F{7], e,X) = \\V- U{v, e, A)f + W.iv) + We{e) + Wx{X), v>0,X>0, 
where the regularizing penalty functionals are of the form 

M M 

(5.24) iy,(r/) = aoJ^Vi + (^iJ2Y. '^^ " ^i'"' 

i=l i=l j&Afi 

M M 

(5.25) Wo{e) = + PiY.J2 1^^'" - ^^'^ I'' 

i=l i=l j&Ni 

(5.26) W^{X) = /?2 (log(A) + u-^ \og{Xf) 

and Mi denotes the usual 4-point nearest neighborhood system for pixel i in the lattice. 

Our objective is to minimize the functional (j5.23p by gradient based optimization meth- 
ods. Here we face the difficulty due to the positivity constraints. To take the positivity 
constraint into account we employ the interior point search method [7j. In the interior point 
search the original constrained problem ()5.23|) is replaced by a sequence of augmented un- 
constrained problems of the form 

(5.27) F, if], 9, A) = F(77, 9, A) + W^^ (r^) 
where w!^\ri) is a penalty functional of the form 

M 

(5.28) w!^\v)=^,Y.- 

i=i 

and {^j} is a sequence of decreasing positive parameters such that ^ as j — > oo. Using 
a suitably chosen sequence of penalty functionals the solutions of the unconstrained 

problems converge (asymptotically) to the solution of the original constrained problem. 
The positivity constraint for A can be taken care with similar techniques. However, it is 
our experience that the positivity constraint was not needed for A. 

For the minimization of the functionals ()5.27|) we employ the Gauss-Newton optimization 
method with an explicit line search algorithm. 
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5.2. Results In this section, we evaluate the performance of the proposed method with 
three different test cases. The first test case is EIT data from an elhpse domain Q, in the 
second test case we consider an elhpse domain with a sharp cut and in the last test case 
the domain is a smooth Fourier domain which has some resemblance with the cross section 
of the human body. In all of these cases, we use the unit disk as the model domain Qm- 

In the simulations, we assume an EIT system with J = 16 electrodes. In each of the test 
cases, the electrodes were located at approximately equally spaced positions at the exterior 
boundary dQ of the target domain Q. The size of the electrodes were chosen such that the 
electrodes covered approximately 50% of the boundary dfl. 

The EIT measurements were simulated using the usual adjacent pair drive data acquisi- 
tion method. In the adjacent drive method, currents +1 and —1 are injected through two 
neighboring electrodes, say electrodes e„ and e„+i, and current through other electrodes is 
zero. The voltages are measured between all J pairs of neighboring electrodes. However, 
three of these measurements are typically neglected since they include either one or both 
of the current feeding electrodes Cn or Cn+i- The rationale behind this is that the electrode 
contact impedances Zj are usually not known accurately. The possible errors in the contact 
impedance values cause a systematic error between the measured voltage and the forward 
model for the measurement made on the current feeding electrodes, and this error causes 
artefacts to the numerical reconstruction, see e.g. [5]. Thus, with the adjacent pair drive 
method each partial measurement consist of L = J — 3 voltage measurements and we have 
yU) g M"^^^. This data acquisition process is then repeated for all the J pairs of adjacent 
electrodes, leading to total of = J(J — 3) voltage measurements for one EIT experiment. 
Thus, with the J = 16 electrode system we have V G 

The simulated EIT measurements were computed using the isotropic EIT model and 
the finite element method. To simulate measurement noise, we added Gaussian random 
noise with standard deviation of 1% of the maximum value of the simulated voltages to 
the data. In all of the following test cases we used value zi = 1 for the electrode contact 
impedances. These were assumed known in the inverse problem. 

The results for the first test case are shown in Figures ^|21 The target conductivity is 
shown in the top left image in Figure ^ The target domain Q is an ellipse with main axes 
1.25 in the horizontal direction and 0.8 in the vertical direction. For the simulation of the 
EIT measurements, the domain was discretized into a finite element mesh that consisted 
of 1256 nodal points and 2350 triangular elements. 

The reconstruction of the conductivity 7 with isotropic EIT model in the correct domain 
Q is shown in the top right image in Figure H The reconstruction was obtained by using 
similar optimization techniques that are explained in the previous section. However, in the 
case of isotropic model the unknown parameter vector is the conductivity vector 7 G M*^ 
and the optimization functionals for the interior point search can be written as 

(5.29) H,{^) = \\V - f/(7)f + W,{^) + W^\^), 

where U{'y) denotes the forward problem for the isotropic model, 1^7(7) and w!^\'y) are 
defined by equations ()5.24|) and ()5.28|) . respectively. To compute the reconstruction in the 
top right image in Figure ^ the domain Q was triangulated to a finite element mesh that 
consisted of 2326 elements with 1244 nodal points. The conductivity was represented in a 
lattice of M = 451 pixels (i.e., 7 G M^^^). The regularization parameters for the penalty 
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Figure 1: Test case with EIT data from an ellipse domain Q. The main axes of the ellipse 
were 1.25 in horizontal direction and 0.8 in the vertical direction. Top left: Simulated 
conductivity distribution 7. Top right: Reconstruction of 7 with isotropic EIT model in 
the correct domain Q. Bottom left: Reconstruction of 7 with isotropic model in incor- 
rectly modeled geometry. The reconstruction domain Q^n was the unit disk. Bottom right: 
Reconstruction of rj with the uniformly anisotropic model in the same unit disk geometry. 

functional W^{'j) in equation ()5.29|) were ao = 10~* and ai = 10~^. When computing 
the reconstruction in the correctly modeled geometry, the interior point search was kept 
inactive (i.e., the sequence of interior point search parameters were all zeros). The 
conductivity vector was initialized to constant value of one in the optimization process. 
The Gauss-Newton optimization algorithm was iterated until convergence was obtained. 

The image in the bottom left in Figure ^ shows the reconstruction of the conductivity 
with the isotropic model in incorrectly modeled geometry Qm- In this case, the computa- 
tional domain Q^n "was the unit disk which was triangulated to 2190 elements with 1176 
nodal points. The conductivity parameters were represented in a lattice of M = 437 pixels 
(i.e., 7 G M^'^^). The regularization parameters for the penalty functional W^{'y) in equation 
(j5.29|) were ao = 10^^ and ai = 2 ■ 10~^. The sequence of interior point search parameters 
{^j} "were from 2 ■ 10""^ to 5 ■ 10""^. The constant vector 7 = 1 G M^^^ was used as the initial 
guess in the Gauss-Newton optimization. 

The image in the bottom right in Figure Q shows the reconstruction of rj with the 
uniformly anisotropic model in incorrectly modeled geometry 0,^- Here, by the solution 
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Figure 2: Test case with EIT data from an ellipse domain ^l. The main axes of the ellipse 
were 1.25 in horizontal direction and 0.8 in the vertical direction. Left: Reconstruction of 
the anisotropy angle parameter 6 in the incorrectly modeled geometry. The computational 
domain Qrn was the unit disk. Right: Evolution of the anisotropy parameter A during the 
Gauss-Newton iteration. 

in the uniformly anisotropic model we mean the optimal solution of the form ()1.5|1 of 
the minimization problem. The reconstruction was obtained by minimizing sequence of 
optimization functionals of the form ()5.27|) . The reconstructed angle parameter 6 is shown 
in the left image in Figure and the evolution of the parameter A during the iteration is 
shown in the right image in Figure |21 The computational domain was the unit disk. 
The finite element triangularization and the amount of the image pixels were the same as 
in the isotropic case in bottom left image in Figure Q Thus, the unknowns in the inverse 
problem are rj G M"^^^, 6 G M^^^ and A G M. The parameters for the regularizing penalty 
functionals Wri{ri) in equation ()5.24p were ao = 10~^ and ai = 10~^. The parameters for the 
penalty functionals We{9) and W\{X) in equations ()5.25H5.25j) were Po = 10~^, Pi = 5 - 10"^ 
and ^2 = 0, respectively. The sequence of interior point search parameters {^j} was from 
1 ■ 10~^ to 1 • 10~^^. The Gauss-Newton optimization was started from the constant values 
?7 = 1 G M^^^, ^ = G M^'^^ and A = 1 which correspond to isotropic unit conductivity. 

The results for the second test case are shown in Figure El The simulated conductivity 
distribution is shown in the top left image. In this case the domain f2 is a truncated ellipse 
with main axes 1.1 in the horizontal direction and 0.9 in the vertical direction, respectively. 
For the simulation of the EIT measurements, the domain was divided to a finite element 
mesh of 2383 triangular elements with 1240 nodes. 

The top right image in Figure El shows the reconstruction of the conductivity with the 
isotropic model in the correct geometry. For the reconstruction, the domain Q was divided 
to a finite element mesh of 2337 triangular elements with 1217 nodes and the conductivity 
was represented on a lattice of M = 455 pixels. Thus, the unknown parameter vector 
was 7 G M^^^. The regularization parameters for the penalty functional W^Ij) in equation 
()5.29|) were ao = 10~^ and ai = 10~^. The sequence of interior point search parameters 
{^j} was all zeros. The Gauss- Newton optimization was started from the constant unit 
conductivity. 

The bottom left image in Figure El shows the reconstructed conductivity with the 
isotropic model in the incorrectly modeled geometry. The reconstruction domain Qrn was 
the unit disk. The finite element mesh and pixel lattice were the same that were used 
for the unit disk in Figured Thus, the unknown conductivity vector was 7 G M'*^^. The 
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Figure 3: Test case with EIT data from a truncated ellipse domain Q with main axes 
a = 1.1 and b = 0.9. Top left: Simulated conductivity distribution 7. Top right: Recon- 
struction of the conductivity 7 with isotropic model in the correct geometry Q. Bottom 
left: Reconstruction of 7 with the isotropic model in incorrectly modeled geometry. The re- 
construction domain Qrn was the unit disk. Bottom right: Reconstruction of the parameter 
1] with the uniformly anisotropic model in the same unit disk geometry. 

parameters for the regularizing penalty functional W.y(7) were ao = 10~^ and ai = 10"'*, 
and the sequence of interior point search parameters {^j} was from 10~^ to 10~^. The 
constant unit conductivity was used as the initial guess in the optimization. 

The bottom right image in Figure El shows the reconstruction of 77 with the uniformly 
anisotropic model in the incorrectly modeled geometry. The computational domain was 
the unit disk with the same discretization that was used in Fig. ^ Thus, the unknowns 
were rj G Mf^^'^, 9 G M^^'' and A G M. The parameters for the regularizing penalty functionals 
Wr^ijj) in equation (j5.24j) were = 10~^ and ai = 10~'^. The parameters for the penalty 
functionals Wg{9) and Wx{X) in equations (I5.25ll5.2(il) were po = 10-^ A = 5 ■ 10"^ and 
P2 = 0, respectively. The sequence of parameters {^j} was from 10~^ to 10~^^. The 
initializations for the parameters in the Gauss-Newton optimization were the constant 
values T] = le M^^^ ^ = G M^^^ and A = 1. 

The results for the last test case are shown in Figure |3] In this case, the target domain Q 
is bounded by a smooth Fourier boundary dfl. The true isotropic conductivity distribution 
within the domain Q is shown in the top left image in Figure El For the simulation of the 
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Figure 4: Test case with EIT data from an arbitrary domain fi. Top left: True conductivity 
distribution 7. Top right: Reconstruction of the conductivity 7 with isotropic model in 
the correct geometry Q. Bottom left: Reconstruction of 7 with the isotropic model in 
incorrectly modeled geometry. The reconstruction domain Qm was the unit disk. Bottom 
right: Reconstruction of the parameter rj with the uniformly anisotropic model in the same 
unit disk geometry. 

EIT measurements, the domain Q was divided to a mesh of 2316 triangular elements with 
1239 nodes. 

The reconstruction of the conductivity 7 with the isotropic model in the correct geometry 
Q is shown in the top right image in Figure EJ The domain was divided to a mesh of 2200 
triangular elements with 1181 nodes for the image reconstruction process. The number of 
pixels was M = 446 for the representation of the conductivity image (i.e. 7 G M"^^^). The 
regularization parameters for the penalty functional W^{'~f) were ao = 10~^ and ai = 10~^ 
and the sequence of parameters {^j} was all zeros. The constant unit conductivity was 
used as the initial guess in the Gauss-Newton optimization algorithm. 

The reconstruction of the conductivity 7 with the isotropic model in the incorrectly 
modeled geometry is shown in the bottom left image in Figure 0] The reconstruction 
domain Qm was the unit disk. The finite element mesh and the pixel lattice were the 
same that were used in Figures Thus, the parameter vector in the inverse problem 
was 7 G W^^"^ . The parameters in the penalty functional W^{'y) were ao = 10~^ and 
ai = 2 ■ 10~^, and the sequence of parameters {^j} was from 2 • 10^^ to 5 • 10~^. The 
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constant unit conductivity was used as the initial guess in the optimization. 

The reconstruction of 7] with the uniformly anisotropic model in the incorrectly modeled 
geometry is shown in the bottom right image in Figure 0] The reconstruction domain Q^n 
was the unit disk with the same discretization as in Figures. Thus, the unknown 

parameter vectors were r] G M^^"^, 6 G M^'^'^ and A G M. The parameters for the regularizing 
penalty functional Wr^{ri) in equation (j5.24|) were = 10^^ and ai = 10~^. The parame- 
ters for the penalty functionals Wo{9) and Wx{\) in equations (|5. 25115. 25|) were Pq = 10~^, 
Pi = 5 ■ 10~^ and P2 = 0, respectively. The sequence of parameters was from 10~^ to 
10~^^. The initializations for the image parameters were the constant values 17 = 1 G R^^^, 
e = Oe M^37 and A = 1. 

6. Discussion As can be seen from Figures ^111 the proposed approach gives good re- 
sults. In all test cases, the traditional reconstructions with the isotropic model are erroneous 
when the imaging geometry is modeled incorrectly. The effects of erroneous geometry are 
seen in the reconstructions as distortions and severe artefacts, especially near the boundary. 
On the other hand, the reconstructions of rj with the uniformly anisotropic model in the 
same erroneous geometry are clear of these artefacts and represent a deformed picture of 
the original isotropic conductivity. These results indicate that the proposed method offers 
an efficient tool to eliminate the difficulties that arise from inaccurately known geometry 
in practical EIT experiments. 
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